#include <stdio.h>
#include <stdlib.h>
#include <getopt.h>
#include <time.h>
#include <sys/time.h>
#include <unistd.h>

#define NUMCOEF 1024 
#define NUMTHREAD 256
#define NUMBLOCK (1024/NUMTHREAD)*8

#define ERROR(fmt, args ...)	\
	printf("[Error] %s:%d\t" fmt, __FUNCTION__, __LINE__, ##args);

static int SIZE;			// Matrix Size
static int ROW_SIZE;
double dL2Norm;
double* gMatrixA;
double* gMatrixB;
double* gCoeffi;

double* dCoeffi;
double* gMatrixdA;

struct option input_option[] = { 
	{"size", 1, 0, 's'},
	{0, 0, 0, 0}
};
struct timeval startTime, endTime;

void show_matrix(double* gMatrix)
{
	int i, j;
	for(i=0; i<1024; ++i){
		printf("%d\t", i);
		for(j=0; j<16; ++j){
			printf("%f\t", gMatrix[i*(SIZE+1)+j]);
		}
		printf("\n");
	}
	printf("\n");
}

int init()
{
	int i, j;

	// Make Input matrix
	gMatrixA = (double*)malloc(sizeof(double)*SIZE*(SIZE+1));
	if(gMatrixA == NULL)	return -1;
	gMatrixB = (double*)malloc(sizeof(double)*SIZE*(SIZE+1));
	if(gMatrixB == NULL)	return -1;

	gCoeffi = (double*)malloc(sizeof(double)*SIZE);

	// For test
	
	srand48(time(0));
	for(i=0; i<SIZE; ++i){
		for(j=0; j<SIZE+1; ++j){
			gMatrixA[i*ROW_SIZE+j] = gMatrixB[i*ROW_SIZE+j] = drand48();
		}
	} 
	return 0;
}

int fini()
{
	free(gMatrixA);
	free(gMatrixB);

	return 0;
}

void set_option(int argc, char** argv)
{
	int opt;
	if(argc!=3 && argc!=5 ){
		printf("Usage: %s [options]\n", argv[0]);
		printf("\t-s\t--size\t\t Size of Matrix\n");
		exit(0);
	}
	do{ 
		opt = getopt_long(argc, argv,"s:t:", input_option, NULL); 
		switch(opt)
		{
			case 's':
				SIZE = atoi(optarg);
				ROW_SIZE = SIZE+1;
				break;
			case -1:
				break;
			default:
				abort();
		}
	}while(opt!=-1);
}

void print_result()
{
	printf("\nSize: %d\n", SIZE);
	printf("L2Norm: %f\n", dL2Norm);
	printf("Time Ellapsed: %f msec\n", (endTime.tv_sec - startTime.tv_sec)*1000 + (endTime.tv_usec - startTime.tv_usec)/1000.0);
}


__global__ void coeffi_calc(double* matrixdA, double* dCoeffi, int ROW_SIZE, int round)
{
	int tx = threadIdx.x;
	int bx = blockIdx.x;
	matrixdA[round * ROW_SIZE + bx * NUMTHREAD + tx] /= matrixdA[round * ROW_SIZE + round];
	dCoeffi[bx * NUMTHREAD + tx] = matrixdA[ROW_SIZE * (bx * NUMTHREAD + tx) + round];
	__syncthreads();
}

__global__ void gaussian_elim(double* matrixdA, double* dCoeffi, int ROW_SIZE, int SIZE, int round)
{
	int i, j;
	int bx = blockIdx.x;
	int tx = threadIdx.x;

	__shared__ double sm_target_row[NUMTHREAD];
	__shared__ double sm_coeffi[NUMCOEF];

	sm_target_row[tx] = matrixdA[round * ROW_SIZE + bx * NUMTHREAD + tx];
	__syncthreads();
	for(i=0;i<SIZE/NUMCOEF;++i) {	
		sm_coeffi[4*tx] = dCoeffi[i * NUMCOEF + 4*tx];
		sm_coeffi[4*tx + 1] = dCoeffi[i * NUMCOEF + 4*tx + 1];
		sm_coeffi[4*tx + 2] = dCoeffi[i * NUMCOEF + 4*tx + 2];
		sm_coeffi[4*tx + 3] = dCoeffi[i * NUMCOEF + 4*tx + 3];
		__syncthreads();
		for(j=0;j<NUMCOEF;++j) { 
			if(i * NUMCOEF + j == round)
				continue;	
			matrixdA[i * NUMCOEF * ROW_SIZE + j * ROW_SIZE + bx * NUMTHREAD + tx] -= sm_coeffi[j]*sm_target_row[tx];
		}
		__syncthreads();
	}
	return;
}

__global__ void normalize(double* matrixdA, int ROW_SIZE)
{
	int bx = blockIdx.x;
	int tx = threadIdx.x;
	int order = (bx * NUMTHREAD) + tx;
	matrixdA[order * ROW_SIZE + ROW_SIZE -1] /= matrixdA[order * ROW_SIZE + order];
	matrixdA[order * ROW_SIZE + order] = 1.0;
}



int device_malloc()
{
	// Copy to device
	if(cudaSuccess != cudaMalloc(&gMatrixdA, sizeof(double)*SIZE*(SIZE+1))){
		ERROR("Fail to cudaMalloc()");
		return -1;
	}
	if(cudaSuccess != cudaMalloc(&dCoeffi, sizeof(double)*SIZE)){
		ERROR("Fail to cudaMalloc()");
		return -1;
	}
	cudaMemcpy(gMatrixdA, gMatrixA, sizeof(double)*SIZE*(SIZE+1), cudaMemcpyHostToDevice);

	return 0;
}


int main(int argc, char** argv)
{
	set_option(argc, argv);

	// Host Initialize
	if(init()<0){
		ERROR("Fail to init()"); 
		return -1;	
	}
    //show_matrix(gMatrixA);
	if(device_malloc() <0){
		ERROR("Fail to device_malloc()");
		return -1;
	}
	dim3 cBlock(NUMTHREAD, 1);
	dim3 cGrid(NUMBLOCK, 1);

	dim3 gBlock(NUMTHREAD, 1);
	dim3 gGrid(NUMBLOCK, 1);

	gettimeofday(&startTime, NULL);
	cudaMemcpy(gMatrixdA, gMatrixA, sizeof(double)*SIZE*(SIZE+1), cudaMemcpyHostToDevice);

	int i,j;
	for(i=0;i<SIZE;++i) {
		coeffi_calc<<<cGrid, cBlock>>>(gMatrixdA, dCoeffi, ROW_SIZE, i);
		gaussian_elim<<<gGrid, gBlock>>>(gMatrixdA, dCoeffi, ROW_SIZE, SIZE, i);
	}
	//normalize<<<gGrid, gBlock>>>(gMatrixdA, ROW_SIZE);

	// Copy to Host
	cudaMemcpy(gMatrixA, gMatrixdA, sizeof(double)*SIZE*(SIZE+1), cudaMemcpyDeviceToHost);
	gettimeofday(&endTime, NULL);

	show_matrix(gMatrixA);

	// Free device memory
	cudaFree(gMatrixdA);

	// show_matrix(gMatrixA);

	print_result();

	// Host Finalize
	fini();
	return 0;
}
